Empirical Generator

In this section, we will see how to construct a generator from data. This is useful when we don't know the underlying generator. There is an inherent limitation in constructing generators that wasn't present when constructing the transfer operator: The data comes at a fixed timescale, thus the generator will only be known up to the timescales present in the timeseries. Nonetheless it is possible to construct a generator from the timeseries under a few assumptions.

For example, suppose that we have the Markov chain with three states

markov_chain = [1 1 2 2 2 1 3 1 1 2 2]

1×11 Matrix{Int64}:
 1  1  2  2  2  1  3  1  1  2  2

and we want to estimate the transfer operator from this data. We can do this by using the function generator which is imported from the MarkovChainHammer.TransitionMatrix module. This function takes a Markov chain as input and returns the empirical transition matrix of the Markov chain. We assume that the time intervals are uniform in time with timestep $dt = 1$

using MarkovChainHammer.TransitionMatrix: generator
dt = 1.0
Q = generator(markov_chain; dt = dt)
Q

3×3 Matrix{Float64}:
 -0.6   0.4   1.0
  0.4  -0.4   0.0
  0.2   0.0  -1.0

This matrix is constructed by counting the amount of time spent in a state to construct the diagonal entries of the matrix, i.e., the holding times. We can manually import the holding times to understand the construction of the $Q$ matrix, from the MarkovChainHammer.TransitionMatrix function holding_times

using MarkovChainHammer.TransitionMatrix: holding_times
state_holding_times = holding_times(markov_chain, 3; dt = dt)

3-element Vector{Vector{Float64}}:
 [2.0, 1.0, 2.0]
 [3.0, 2.0]
 [1.0]

The first argument to the function is the markovchain, the second argument is the total number of states, and the keyword argument is the time step size associated with the entries of the markovchain.

We see, from the markov chain, that state 1 spends 2 time units in state 1, followed by 1 time unit, then follow by 2 time units. State 2 spends 3 time units repeating itself, then two time units. And lastly, state 3 is only observed once for one time unit.

The diagonals of the $Q$ matrix are given by taking the average of the holding times, then taking the reciprocal. We can verify this manually by importing the mean function from the Statistics package.

using Statistics
1 ./ mean.(state_holding_times)

3-element Vector{Float64}:
 0.6
 0.4
 1.0

We can verify the off diagonal terms as well. We just need to track the number of times that a given state was exited. For example, when leaving state 1 we observe $1 \rightarrow 2$, read as 1 transitions to 2, $1 \rightarrow 3$, and $1 \rightarrow 2$. Thus going from state state 1 to state 2 is twice as likely as going from state 1 to state 3 and the associated probabilities are $2/3$ for transitioning from 1 to 2 and $1/3$ for transitioning from 2 to 3. These two probabilities then get multiplied by the reciprocal empircal holding time for being in state 1, which in this case is $0.6$ to yield the first column of the $Q$ matrix

Q[:, 1]

3-element Vector{Float64}:
 -0.6
  0.39999999999999997
  0.19999999999999998

For the other two states we observe that we only see state 2 transitioning to state 1 and state 3 only ever transitions to state 1, thus the matrix entries are $Q_{32} = Q_{23} = 0$. The requirement that the columns must sum to zero then determines the other entries (or noticing that the probability of going to state 1 after having left state 2 is 1 and similarly for state 3).

This brings us to our first warning about constructing the empirical generator. It is most meaningful when there is enough timeseries data to stay within a given state for more than one "timestep".

Our last comment is that changing $dt$ amounts to rescaling the generator. That is to say,

Q1 = generator(markov_chain; dt = 1.0)
Q2  = generator(markov_chain; dt = 2.0)
Q1 - 2.0 * Q2

3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

As a last comment we mention that one can also construct the transfer operator, then take the matrix logarithm, and then divide by $dt$ to get another estimate of the generator; however, the resulting matrix no longer has an "infinitesimal" probabilistic interpretation, as the following example shows

using MarkovChainHammer.TransitionMatrix: perron_frobenius
dt = 1.0
ℳ = perron_frobenius(markov_chain)
log(ℳ) / dt

3×3 Matrix{ComplexF64}:
 -0.612402+1.03113im    0.299568-0.241032im    0.665466-3.22739im
  0.479309-0.385651im  -0.335089+0.0901481im   0.284164+1.20707im
  0.133093-0.645477im  0.0355205+0.150884im   -0.949629+2.02032im

The columns still sum to zero, but the entries of the matrix no longer have an interpretation as holding times or probabilities. This example also shows why generators have a much more limited and special structure as compared to the transfer operator. The requirement that one generates probabilities over infinitesimal steps imposes severe restrictions. We may lift these restrictions if we are willing to give up computing an infinitesimal generator and live with generators that are applicable only over a finite timescale. This fact reflects a well-known theorem in numerical analysis that positivity preserving linear operators are inherently lower-order.

An alternative way to construct the generator from the perron-frobenius operator is to take the difference between the identity matrix and the perron-frobenius operator and then divide by dt

using LinearAlgebra
dt = 1.0
Q = (ℳ - I) / dt

3×3 Matrix{Float64}:
 -0.6   0.25   1.0
  0.4  -0.25   0.0
  0.2   0.0   -1.0

where we had to use the LinearAlgebra package to import the identity matrix. This calculation makes use of the relation that the generator is the infinitesimal generator of the perron-frobenius operator

\[\mathcal{M} = e^{dt Q} \approx \mathbb{I} + dt Q \]

This is only good if the timesteps are sufficiently small.